#-*- coding:utf-8 -*-
import time, math, random
import numpy as np

# python标准库里的优先队列
from multiprocessing import Queue
from queue import PriorityQueue
# python标准库里的队列
from collections import deque

from igraph import *
# from prettytable import PrettyTable

#表示无穷大
DBL_MAX = 1.0e6
#放大因子
LARGE_COEFF = 1e3
#增加虚拟源点以及虚拟的出边
def add_virtual_source(dg):
    #查找所有入度等于0的节点
    nodes = dg.vs.select(_indegree=0)
    if len(nodes) == 0:return -1
    #增加一个虚拟源点
    dg.add_vertex()
    sn = dg.vs.indices[-1]
    #增加虚拟分支(虚拟源点-->)
    for u in nodes:
        dg.add_edge(sn, u.index)
    #返回新增的虚拟源点
    return sn
#增加虚拟汇点以及虚拟的入边
def add_virtual_target(dg):
    #查找所有出度等于0的节点
    nodes = dg.vs.select(_outdegree=0)
    if len(nodes) == 0:return -1
    #增加一个虚拟汇点
    dg.add_vertex()
    tn = dg.vs.indices[-1]
    #增加虚拟分支(-->虚拟汇点)
    for u in nodes:
        dg.add_edge(u.index, tn)
    #返回新增的虚拟汇点
    return tn
#e是一个igraph.Edge类型变量
#计算风机特性曲线的多项式值
def h0(e):
    q, a0, a1, a2 = e['q'],e['a0'],e['a1'],e['a2']
    #计算多项式
    #多项式函数用法参考：numpy-html-1.7.0/reference/generated/numpy.polyval.html#numpy.polyval
    return np.polyval([a2, a1, a0], q)
#e是一个igraph.Edge类型变量
#计算风机特性曲线的1阶导
def h1(e):
    q, a0, a1, a2 = e['q'],e['a0'],e['a1'],e['a2']
    #计算多项式的一阶导
    #多项式函数用法参考：numpy-html-1.7.0/reference/generated/numpy.polyval.html#numpy.polyval
    return np.polyval([2*a2, a1], q)
#e是一个igraph.Edge类型变量
#分支的摩擦风阻
def R0(e):
    return e['r']
#e是一个igraph.Edge类型变量
#分支的总风阻(包括调节风阻)
def R(e):
    return e['r']+e['adjust_r']
#e是一个igraph.Edge类型变量
#分支的调节风阻(通过调节压力，也就是压差来计算)
def dr(e):
    return 0.0 if abs(e['q']) < 1e-4 else e['adjust_h']/(e['q']**2)
#e是一个igraph.Edge类型变量
#分支的阻力h=r*q^2
def f0(e):
    # print(e['id'],R(e),e['q'])
    return R(e)*e['q']*abs(e['q'])
#e是一个igraph.Edge类型变量
#分支阻力的1阶导
def f1(e):
    return 2*abs(e['q']*R(e))
#e是一个igraph.Edge类型变量
#是否虚拟分支
def is_zero_edge(e):
    return abs(100*e['r']) < 1e-6
#e是一个igraph.Edge类型变量
#是否固定风量分支
def is_fix_edge(e):
    return abs(e['fq']) > 0
#e是一个igraph.Edge类型变量
#是否风机分支(分支上有风机)
def is_fan_edge(e):
    return not (e['a0']==0 and e['a1']==0 and e['a2']==0)
#创建图的拓扑关系(节点和分支,以及实际id与igraph内部编号之间的映射)
def build_graph(edges, dg):
    #收集所有的节点
    sIDs = [str(data['s']) for data in edges]
    tIDs = [str(data['t']) for data in edges]
    sIDs.extend(tIDs)
    #构造节点集合(利用set去除重复编号,并排序,然后再转换成list)
    nIDs = set(sIDs)
    #构造节点编号到igraph内部编号的映射关系
    nodes_map = dict(zip(nIDs, range(len(nIDs))))
    #有向图增加节点(igraph节点内部编号从0开始)
    dg.add_vertices(len(nodes_map))
        
    #构造分支编号到igraph内部编号的映射关系
    edges_map = {}
    #分支的实际编号
    e_id = 0
    for data in edges:
        eID, u, v = str(data['id']), str(data['s']), str(data['t'])
        u, v = nodes_map[u], nodes_map[v] #转换成igraph内部的编号
        # 不考虑重复边!!!
        if eID not in edges_map:
            dg.add_edge(u, v)  #增加新分支(igraph分支内部编号从0开始)
            edges_map[eID] = e_id #记录分支编号到igraph内部编号的映射关系
            e_id = e_id + 1
    return nodes_map, edges_map
def set_graph_ids(dg, nodes_map, edges_map):
    #设置igraph节点的id属性数据
    for u in nodes_map:
        i = nodes_map[u]
        dg.vs[i]['id'] = u
    #设置igraph分支的id属性数据
    for e in edges_map:
        i = edges_map[e]
        dg.es[i]['id'] = e
#初始化图的属性默认值
def init_graph_property(dg):
    dg.es['id']=str(0)
    dg.es['a0']=0.0
    dg.es['a1']=0.0
    dg.es['a2']=0.0
    # 风机多个解时，用于选择斜率小于0的值调节
    dg.es['fan_c']=200.0
    dg.es['r']=0.0
    dg.es['r0']=0.0
    dg.es['q']=0.0
    dg.es['fq']=0.0
    dg.es['adjust_r']=0.0
    dg.es['adjust_h']=0.0
    dg.es['density']=1.2
    #节点数据
    dg.vs['id']=str(0)
    dg.vs['p']=0.0
    dg.vs['x']=0.0
    dg.vs['y']=0.0
    dg.vs['z']=0.0
    #图数据
    dg['Q'] = 100.0
    dg['q_precise'] = 1e-4 # 修正风量精度
    dg['h_precise'] = 1.0 # 回路阻力代数和精度
    dg['maxCount'] = 1000
#如果属性值是None,则给它设定默认值
def make_defalut_value(obj, property_name, property_value):
    if obj is None or property_name not in obj.attributes() or property_value is None:
        return
    elif obj[property_name] is None:
        obj[property_name] = property_value
#遍历图中的所有分支和节点,如果它的属性值为None,则给默认值
def default_graph_property(dg):
    for e in dg.es:
        make_defalut_value(e, 'id', str(0))
        make_defalut_value(e, 'a0', 0.0)
        make_defalut_value(e, 'a1', 0.0)
        make_defalut_value(e, 'a2', 0.0)
        make_defalut_value(e, 'fan_c', 200.0)
        make_defalut_value(e, 'r', 0.0)
        make_defalut_value(e, 'r0', 0.0)
        make_defalut_value(e, 'q', 0.0)
        make_defalut_value(e, 'fq', 0.0)
        make_defalut_value(e, 'adjust_r', 0.0)
        make_defalut_value(e, 'adjust_h', 0.0)
        make_defalut_value(e, 'density', 1.2)
    for u in dg.vs:
        make_defalut_value(u, 'id', str(0))
        make_defalut_value(u, 'p', 0.0)
        make_defalut_value(u, 'x', 0.0)
        make_defalut_value(u, 'y', 0.0)
        make_defalut_value(u, 'z', 0.0)
    make_defalut_value(dg, 'Q', 100.0)
    make_defalut_value(dg, 'q_precise', 1e-4)
    make_defalut_value(dg, 'h_precise', 1.0)
    make_defalut_value(dg, 'maxCount', 1000)
#构建通风网络
def build_network(graph_datas, dg):
    if 'edges' not in graph_datas:
        return False
    elif len(graph_datas['edges']) == 0:
        return False
    #从词典中提取分支数据
    edges = graph_datas['edges']
    #从词典中提取节点数据
    nodes = graph_datas['nodes']
    #构建有向图并返回节点编号、分支编号与igraph内部编号的映射关系
    nodes_map, edges_map = build_graph(edges, dg)
    #初始化所有的分支数据默认值
    #如果在dg.add_vertices()或dg.add_edge()之前设置的属性会出现一个问题:
    #    新增的分支的属性值为None
    #所以这里为了保险期间,在分支和节点都被add完事之后,再给节点和分支挂载属性
    init_graph_property(dg)
    #设置分支和节点的id属性
    set_graph_ids(dg, nodes_map, edges_map)
    #设置分支预定义数据
    for data in edges:
        i = edges_map[str(data['id'])]  #igraph分支的内部编号
        dg.es[i]['r'] = data.get('r', 0.0)  #设置分支的风阻(预定义)
        dg.es[i]['r0'] = data.get('r0', 0.0)  #设置分支的初始风阻(预定义)
        dg.es[i]['q'] = data.get('q', 0.0)  #设置分支的风量(预定义)
        dg.es[i]['adjust_r'] = data.get('adjust_r', 0.0)  #设置分支的风量(预定义)
        # 判断风阻是否为负数，风阻为负数则不需要进行下一步的解算等操作
        # by huangde 2020/1/6
        if dg.es[i]['adjust_r'] + dg.es[i]['r'] < 0:
            print('分支[%s]风阻为负值'%data['id'])
            return False
    # 设置节点预定义数据
    for data in nodes:
        i = nodes_map[str(data['id'])]  #igraph节点的内部编号
        dg.vs[i]['p'] = data.get('p', 0.0)  #设置节点的压能(预定义)
        dg.vs[i]['p_used'] = data.get('p_used', False)  #是否使用预定义的节点压能
        # 节点压能计算是否加上风机的负压，全局变量
        if 'p_used_fans' in data:
            dg.vs[i]['p_used_fans'] = data.get('p_used_fans',True)
        else:
            dg.vs[i]['p_used_fans'] = graph_datas.get('p_used_fans',True)

    #设置分支其它数据(权重、费用等)
    if 'edge_weights' in graph_datas:
        for data in graph_datas['edge_weights']:
            i = edges_map[str(data['id'])]  #igraph分支的内部编号
            for name in data:
                if name == 'id':continue
                dg.es[i][name] = data[name]
    #设置分支其它数据(权重、费用等)
    if 'node_weights' in graph_datas:
        for data in graph_datas['node_weights']:
            i = nodes_map[str(data['id'])]  #igraph节点的内部编号
            for name in data:
                if name == 'id':continue
                dg.vs[i][name] = data[name]
    #设置风机数据
    if 'fans' in graph_datas:
        fans = graph_datas['fans']
        for data in fans:
            if data['fan_flag'] == 0:
                continue
            i = edges_map[str(data['e_id'])]  #igraph分支的内部编号
            dg.es[i]['a0'] = data['fan_flag']*data.get('a0', 0.0)
            dg.es[i]['a1'] = data.get('a1', 0.0)
            dg.es[i]['a2'] = data['fan_flag']*data.get('a2', 0.0)
            dg.es[i]['fan_c'] = data['fan_flag']*data.get('fan_c',200.0)
    #设置固定风量数据
    if 'qFixs' in graph_datas:
        qFixs = graph_datas['qFixs']
        for data in qFixs:
            i = edges_map[str(data['e_id'])]  #igraph分支的内部编号
            dg.es[i]['fq'] = data.get('fq', 0.0)
    #设置构筑物
    #构筑物风阻是直接加到分支的调节风阻上的!
    #在vno中构筑物数据并未参与计算过程!
    #todo[by dlj on 2015/2/11]:构筑物的设计后续可能需要重新考虑下!!!
    if 'gates' in graph_datas:
        gates = graph_datas['gates']
        for data in gates:
            i = edges_map[str(data['e_id'])]  #igraph分支的内部编号
            dg.es[i]['adjust_r'] = dg.es[i]['adjust_r'] + data.get('r', 0.0)
    #设置调节参数
    if 'Q' in graph_datas:
        dg['Q'] = graph_datas['Q']
    if 'q_precise' in graph_datas:
        dg['q_precise'] = graph_datas['q_precise']
    if 'h_precise' in graph_datas:
        dg['h_precise'] = graph_datas['h_precise']
    if 'maxCount' in graph_datas:
        dg['maxCount'] = graph_datas['maxCount']
    return True
#独立回路
class IndependentCycle:
    def __init__(self, dg, baseEdge, path):
        self.dg = dg
        self.baseEdge = baseEdge
        self.path = path[:]
        # print('基准分支:e%d, 路径:%s' % (self.baseEdge, str(self.path)))
        self.coef = []
        #计算方向系数
        self.__initDirection()
        # print('方向系数:,self.coef)
        # 判断回路是否可以参与迭代
        self.bCanDoIterate = self.__canDoIterate()
    # 回路的基准分支
    def getBaseEdge(self):
        return self.dg.es[self.baseEdge]
    def getPath(self):
        edges = []
        for i in self.path:
            edges.append(self.dg.es[i])
        return edges
    # 回路中第i条分支的方向(回路按余枝的正方向排列)
    def getCycleCoeff(self, i):
        if i < 0:
            return 0
        elif i == self.baseEdge:
            return 1
        elif i not in self.path:
            return 0
        else:
            return self.coef[self.path.index(i)]
    # 回路阻力代数和(压差)
    def pressureDrop(self):
        return self.__f0() - self.__h0()
    # 回路迭代
    def iterate(self, bCrackFan = True):
        dq, df = 0, 0
        if not self.bCanDoIterate:
            return False, dq, df
        else:
            # 计算回路风量修正量
            dq = self.__deltaQ()
            if abs(dq) > 1e4:
                dq = 10 if dq>0 else -10
            # 特殊处理风机分支
            # 保证风机分支的风量始终在抛物线的右侧(斜率小于0的一侧)
            if bCrackFan:
                dq = self.__crackFan(dq)
            # 修正回路风量
            self.__increaseQ(dq)
            # 回路阻力代数和
            df = self.__f0() - self.__h0()
            return True, dq, df
    def __initDirection(self):
        if len(self.path) == 0:return
        u = self.dg.es[self.baseEdge].target
        d = 1
        for i in self.path:
            s, t = self.dg.es[i].tuple
            if d > 0:
                if s != u:
                    d = -1*d
                    u = s
                else:
                    u = t
            else:
                if t != u:
                    d = -1*d
                    u = t
                else:
                    u = s
            self.coef.append(d)
    def __canDoIterate(self):
        #1) 基准分支是固定风量分支,不参与迭代计算
        #2) 回路中包含固定风量分支,不参与迭代计算
        if self.__isFixQCycle() or self.__hasFixQ():
            return False
        #3) 如果基准分支是大气分支, 且不包含附加阻力, 则不参与迭代计算
        elif self.__isAirCycle() and not self.__hasExtraF():
            return False
        #4) 如果基准分支是未使用的风机分支,且回路中没有附加阻力,则不参与迭代计算
        #elif not self.__hasExtraF():
        #   return False
        return True
    def __isAirCycle(self):
        # 假设1:风阻为0的基准分支是大气分支,即为大气回路
        # 假设2:虚拟汇点和原始网络汇点之间的虚拟分支不可能是基准分支!!!
        return is_zero_edge(self.dg.es[self.baseEdge])
    def __isFixQCycle(self):
        return is_fix_edge(self.dg.es[self.baseEdge])
    #检查回路中是否存在附加阻力(风机)
    def __hasExtraF(self):
        ret = False
        for i in self.path:
            #分支上有附加阻力(风机)
            if is_fan_edge(self.dg.es[i]):
                ret = True
                break
            else:
                #其它的情况,例如自然风压,火风压等等
                #目前还没有考虑!!!
                pass
        return ret
    def __hasFixQ(self):
        ret = False
        for i in self.path:
            if is_fix_edge(self.dg.es[i]):
                ret = True
                break
        return ret
    def __deltaQ(self):
        # 风机动力方向与方向系数c有关
        # sq - c*hf
        # sq表示回路分支阻力代数和,hf表示风机动力(根据特性曲线计算出来的)
        # 而风机分支始终是余支,也即c=1
        F0 = self.__f0()
        H0 = self.__h0()
        F1 = self.__f1()
        H1 = self.__h1()
        # print(-1.0*(F0-H0)/(F1-H1+1e-4))
        return -1.0*(F0-H0)/(F1-H1+1e-4)
    def __increaseQ(self, dq):
        # 修正基准分支风量
        self.__increaseEdgeQ(self.baseEdge, dq)
        # 回路中的其它分支修正风量
        for i, c in zip(self.path, self.coef):
            # e = self.dg.es[i]
            # print 'q%s = q%s+%d*%.3f' % (e['id'],e['id'],c,dq)
            # 修正回路分支风量
            self.__increaseEdgeQ(i, c*dq)
    def __increaseEdgeQ(self, i, dq):
        self.dg.es[i]['q'] = self.dg.es[i]['q'] + dq
    def __f0(self):
        return f0(self.dg.es[self.baseEdge]) + np.sum([f0(self.dg.es[i])*c for i,c in zip(self.path, self.coef)])
    def __f1(self):
        return f1(self.dg.es[self.baseEdge]) + np.sum([f1(self.dg.es[i])*abs(c) for i,c in zip(self.path, self.coef)])
    def __h0(self):
        return h0(self.dg.es[self.baseEdge]) + np.sum([h0(self.dg.es[i])*c for i,c in zip(self.path, self.coef)])
    def __h1(self):
        return h1(self.dg.es[self.baseEdge]) + np.sum([h1(self.dg.es[i])*c for i,c in zip(self.path, self.coef)])
    def __crackFan(self, dq):
        #如果基准分支上有风机,且斜率>0
        # 2019年12月26日，之前给定的修正值200，用于选择右侧值，把这个值写到风机属性中？
        if h1(self.dg.es[self.baseEdge]) > 0:
            # print('基准分支风机斜率大于0,c=',self.dg.es[self.baseEdge]['fan_c'])
            # print(self.dg.es[self.baseEdge]['q'])
            dq = self.dg.es[self.baseEdge]['fan_c']
            return dq
        for i,c in zip(self.path, self.coef):
            if h1(self.dg.es[i]) > 0:
                # print(self.dg.es[i]['id'],'风机斜率大于0,c=',self.dg.es[i]['fan_c'])
                dq = dq + c*self.dg.es[i]['fan_c']
                break
        return dq
#通风网络
class VentNetwork:
    def __init__(self):
        #定义有向图
        self.dg = Graph(directed=True)
        #虚拟源点和汇点以及虚拟大气分支
        self.sn = -1
        self.tn = -1
        self.airEdge = -1
    #返回图对象
    def graph(self):
        return self.dg
    #返回虚拟源点
    def vSource(self):
        return self.sn
    #返回虚拟汇点
    def vTarget(self):
        return self.tn
    # 返回虚拟源汇点
    def vST(self):
        return self.sn, self.tn
    #返回虚拟大气分支
    def vAir(self):
        return self.airEdge
    #增加虚拟的源点和汇点,并附加相关的属性数据
    def addVirtualST(self):
        if self.sn > -1 or self.tn > -1:return False
        #增加虚拟源点和汇点
        self.sn = add_virtual_source(self.dg)
        self.tn = add_virtual_target(self.dg)
        if self.sn < 0 or self.tn < 0:
            return False
        else:
            #初始化虚拟分支和节点的数据
            #这里为了省事,直接遍历所有分支和节点(不单独处理虚拟节点和虚拟分支了)
            default_graph_property(self.dg)
            #虚拟源汇赋予id, 便于后续处理!
            self.dg.vs[self.sn]['id'] = 'S'
            self.dg.vs[self.tn]['id'] = 'T'
            # 特殊处理虚拟分支风量
            for i in self.dg.incident(self.sn, mode=OUT):
                e = self.dg.es[i]
                e['q'] = self.outArcFlow(e.target)
            for i in self.dg.incident(self.tn, mode=IN):
                e = self.dg.es[i]
                e['q'] = self.inArcFlow(e.source)
            return True
    def delVirtualST(self):
        if self.sn < 0 or self.tn < 0:return
        #删除虚拟源点和汇点以及虚拟分支
        self.dg.delete_vertices([self.sn, self.tn])
        self.sn, self.tn = -1, -1
    #增加虚拟的大气分支(虚拟的汇点->虚拟的源点)
    def addAirEdge(self):
        if self.sn < 0 or self.tn < 0 or self.airEdge >= 0:return
        #增加虚拟大气分支
        self.dg.add_edge(self.tn, self.sn)
        self.airEdge = self.dg.es.indices[-1]
        #初始化虚拟分支的数据
        #这里为了省事,直接遍历所有分支和节点(不单独处理虚拟节点和虚拟分支了)
        default_graph_property(self.dg)
        #大气分支的风量等于总风量(源点的出边分支总风量)
        self.dg.es[self.airEdge]['q'] = self.outArcFlow(self.sn)
    #删除虚拟的大气分支(虚拟的汇点->虚拟的源点)
    def delAirEdge(self):
        if self.airEdge < 0:return
        self.dg.delete_edges(self.airEdge)
        self.airEdge = -1
    # 原始网络的源点节点
    def sourceNodes(self):
        # 没有增加虚拟源汇
        if self.sn < 0:
            return self.dg.vs.select(_indegree=0)
        else:
            return [self.dg.es[i].target for i in self.dg.incident(self.sn, mode=OUT)]
    # 原始网络的汇点节点
    def targetNodes(self):
        if self.tn < 0:
            return self.dg.vs.select(_outdegree=0)
        else:
            return [self.dg.es[i].source for i in self.dg.incident(self.sn, mode=IN)]
    def hPath(self, p):
        return np.sum([f0(e) for e in self.dg.es.select(p)])
    def filterVirutalEdges(self, p):
        return [e.index for e in self.dg.es.select(p) if not is_zero_edge(e)]
    def outArcFlow(self, u):
        return np.sum([self.dg.es[i]['q'] for i in self.dg.incident(u, mode=OUT)])
    def inArcFlow(self, u):
        return np.sum([self.dg.es[i]['q'] for i in self.dg.incident(u, mode=IN)])
    def minMaxQ(self):
        edges = [e for e in self.dg.es if not is_zero_edge(e)]
        return min(edges, key=lambda e: abs(e['q'])), max(edges, key=lambda e: abs(e['q']))
    def totalFlow(self):
        if self.sn < 0:
            # 计算多个源点的出边风量之和
            return np.sum([self.outArcFlow(u) for u in self.dg.vs.select(_indegree=0)])
        else:
            # 计算虚拟源点的所有出边风量(a little trick)
            # return self.outArcFlow(self.sn)
            return np.sum([self.outArcFlow(self.dg.es[i].target) for i in self.dg.incident(self.sn, mode=OUT)])
    def totalH(self):
        return self.hPath(mrp(self, self.sn, self.tn, True))
    def totalR(self):
        H = self.totalH()
        Q = self.totalFlow()
        return H/Q**2
    # 是否连通图
    def isConnected(self):
        return self.dg.is_connected(mode=WEAK)
    # 是否有向无环图
    def isDag(self):
        return self.dg.is_dag()
    # 是否有单向回路
    def hasCycles(self):
        # return self.dg.is_connected(mode=STRONG)
        return len(self.SCC()) > 0
    # 查找连通块(弱联通块)
    def CC(self):
        cl = self.dg.components(mode=WEAK)
        comps = []
        for g in cl.subgraphs():
            # 子图的index与原图已经不一致了!
            # 这里只能使用id属性进行返回!
            comps.append([e['id'] for e in g.es])
        # 释放内存
        del cl
        return comps
    # 查找单向回路(强连通块)
    def SCC(self):
        cl = self.dg.components(mode=STRONG)
        # print list(cl)
        comps = []
        for g in cl.subgraphs():
            # 至少2条边才能构成单向回路
            if len(g.es) < 2:continue
            # 子图的index与原图已经不一致了!
            # 这里只能使用id属性进行返回!
            for e in g.es:
                node_path = dfs_path(g,e.target,e.source)
                edge_path = node_path_to_edge_path(g,node_path)
                edge_path.append(e.index)
                edge_path.sort()
                comp = [g.es[i]['id'] for i in edge_path]
                if comp not in comps:
                    comps.append(comp)
        # 释放内存
        del cl
        return comps
    # 统计固定风量的个数
    def countFixQ(self):
        return len([e.index for e in self.dg.es if is_fix_edge(e)])
    # 进风井
    def sourceEdges(self):
        edges = []
        for u in self.dg.vs.select(_indegree=0):
            edges.extend([self.dg.es[i] for i in self.dg.incident(u, mode=OUT)])
        return edges
    # 回风井
    def sinkEdges(self):
        edges = []
        for u in self.dg.vs.select(_outdegree=0):
            edges.extend([self.dg.es[i] for i in self.dg.incident(u, mode=IN)])
        return edges
    # 固定风量分支
    def fixQEdges(self):
        return [e.index for e in self.dg.es if is_fix_edge(e)]
    # 风机分支
    def fanEdges(self):
        return [e.index for e in self.dg.es if is_fan_edge(e)]
    def negativeEdges(self):
        return [e.index for e in self.dg.es if e['q'] < 0]
    # 根据分支id查找分支的内部编号
    def getEdgeById(self, edge_id):
        edge = -1
        for e in self.dg.es:
            if e['id'] == str(edge_id):
                edge = e.index
                break
        return edge
    # 根据分支id查找分支的内部编号
    def getNodeById(self, node_id):
        node = -1
        for v in self.dg.vs:
            if v['id'] == str(node_id):
                node = v.index
                break
        return node
    # 根据风机分支id查找风机分支的内部编号
    def getFanEdgeById(self, edge_id):        
        fan_edges = self.fanEdges()
        if len(fan_edges) == 0:
            return -1
        edge = -1
        # print fan_edges
        for i in fan_edges:
            if self.dg.es[i]['id'] == str(edge_id):
                edge = i
                break
        return edge
    # 调整负风量分支及其方向
    def adjustNegativeEdges(self):
        ne = []
        for e in self.dg.es:
            if e['q'] < 0:
                e['q'] = -1.0*e['q']
                #记录要删除的分支id
                ne.append(e.index)
                #增加新的正方向分支(将原来负方向的分支属性完全复制过去)
                self.dg.add_edge(e.target, e.source, **e.attributes())
        #删除原来的负方向分支
        self.dg.delete_edges(ne)
    # 计算节点压力
    def caclNodePressures(self):
        # 标记节点压力是否已被计算过了
        self.dg.vs['hasPressCacl'] = False
        # 源点的压力记为0,并标记为已计算
        self.dg.vs[self.sn]['hasPressCacl'] = True
        if not self.dg.vs[self.sn]['p_used']:
            self.dg.vs[self.sn]['p'] = 0.0
        # 标记虚拟源点的出边分支的末节点都已经计算过了
        # 并将找到的节点依次添加到队列中,按出队的顺序计算节点压力
        cNodes = deque()
        for i in self.dg.incident(self.sn, mode=OUT):
            e = self.dg.es[i]
            u = e.target # u是一个整数,不是节点对象
            self.dg.vs[u]['hasPressCacl'] = True
            cNodes.append(u)
        # 顺序出队,依次计算节点压力
        while len(cNodes) != 0:
            u = cNodes[0]
            # print u, '出队'
            for i in self.dg.incident(u, mode=OUT):
                e = self.dg.es[i]
                v = e.target # u是一个整数,不是节点对象
                if v == self.sn or v == self.tn or self.dg.vs[v]['hasPressCacl']:
                    continue
                # print v, '进队'
                # 计算分支阻力 和 风机动力
                # print e['id'], f0(e), h0(e)
                # 若使用预定义的节点压能，则不计算
                if self.dg.vs[v]['p_used']:
                    self.dg.vs[v]['p'] = self.dg.vs[v]['p']
                    # print(self.dg.vs[v]['id'], f0(e), h0(e))
                # 根据是否使用风机的负压分别计算
                elif self.dg.vs[v]['p_used_fans']:
                    self.dg.vs[v]['p'] = self.dg.vs[u]['p'] - f0(e) + h0(e)
                else:
                    self.dg.vs[v]['p'] = self.dg.vs[u]['p'] - f0(e)
                cNodes.append(v)
                self.dg.vs[v]['hasPressCacl'] = True
            cNodes.popleft()

# 打印通风网络风量分配结果
def internal_print_network(vnet, msg):
    dg = vnet.graph()
    # print '---------- igraph内部详细信息 ---------'
    # print dg
    print('---------- %s ---------' % msg)
    for e in dg.es:
        # print 'r(e%s)=%.3f\tq(e%s)=%.3f\th(e%s)=%.3f' % (e['id'],e['r'],e['id'],e['q'],e['id'],f0(e))
        # print dir(e)
        print('r(e%d)=%.3f\tq(e%d)=%.3f\th(e%d)=%.3f' % (e.index,e['r'],e.index,e['q'],e.index,f0(e)))
#构造边的权重--用于有上下界的最大流算法
def build_fixq_weight(vnet):
    #计算分支权重
    dg = vnet.graph()
    # 统计固定风量个数
    fixq_count = vnet.countFixQ()
    # print '固定风量个数:',fixq_count
    # 此时通风网络中没有固定风量
    # 如果没有特殊处理,InitFixQ返回0值,风量分配失败!
    # 特殊处理方法:
        # 虚拟源点的第一条出边的风量人为固定一个数值
        # 这就使得InitFixQ方法可以正确执行了
        # 再固定分配完毕后将虚拟源点的第一条出边的固定风量恢复
    if fixq_count == 0:
        # 人为固定第一条出边的风量为初始总风量
        i = dg.incident(vnet.vSource(), mode=OUT)[0]
        dg.es[i]['fq'] = dg['Q']
    #依次给图中的其他分支分配权重
    for e in dg.es:
        e['weight'] = e['fq'] if is_fix_edge(e) else 0.0
    if fixq_count == 0:
        # 恢复第一条出边的固定风量为0
        i = dg.incident(vnet.vSource(), mode=OUT)[0]
        dg.es[i]['fq'] = 0
# 有上下界的最大流算法
def bound_max_flow(dg, s, t, Q, fixTotalQ=False):
    # 在源点s之前再添加一条分支,用来限制总流量Q
    dg.add_vertex()
    s0 = dg.vs.indices[-1]
    dg.add_edge(s0, s)
    e0 = dg.es.indices[-1]
    dg.es[e0]['low'] = Q*1.0 if fixTotalQ else 0.0
    # dg.es[e0]['high'] = Q*1.0
    dg.es[e0]['high'] = Q*1.0 if fixTotalQ else DBL_MAX # UPDATE: 2016/12/24
    # 添加一条t->s的分支,下界0,上界无穷大
    dg.add_edge(t, s0)
    ts = dg.es.indices[-1]
    dg.es[ts]['low'] = 0.0
    dg.es[ts]['high'] = DBL_MAX
    # 分支容量(默认无穷大)
    dg.es['capacity'] = DBL_MAX
    # 上界分支拆分后的分支映射记
    dg.es['lu'] = -1
    # 添加超级源汇
    dg.add_vertex()
    ss = dg.vs.indices[-1]
    dg.add_vertex()
    tt = dg.vs.indices[-1]
    # 注:后续新增的节点和分支属性map不会自动更新,需要手动设置
    # 例如cap[e0]、cap[ts]、arc_filter、lu等
    # 为了减少重复代码,将属性初始化声明转移到添加节点和分支之后
    # 拆分有上下界的分支,其它分支上界默认无穷大
    for e in dg.es:
        if e['low'] > 0 and e['high'] >= e['low']:
            # print 'e%d --> low:%.1f high:%.1f' % (e.index, e['low'], e['high'])
            u, v = e.tuple
            dg.add_edge(ss, v)
            e1 = dg.es[dg.es.indices[-1]]
            dg.add_edge(u, tt)
            e2 = dg.es[dg.es.indices[-1]]
            e['capacity'] = e['high'] - e['low']
            e1['capacity'], e2['capacity'] = e['low'], e['low']
            e1['low'], e2['low'] = 0.0, 0.0
            e1['high'], e2['high'] = DBL_MAX, DBL_MAX
            e1['lu'], e2['lu'] = e.index, e.index
    # for e in dg.es:
    #   print 'e%d --> lu:%d' % (e.index, e['lu'])
    # 执行一次最大流计算
    flow = dg.maxflow(ss, tt, "capacity")
    # 最大流值
    max_flow = flow.value
    # 计算得到的分支流量
    flow_map = flow.flow
    # 判断是否可行流(即额外添加的分支满流,满足下界low)
    # 2种方法:
        # 1) 所有的流量的下界之和是否等于最大流
        # 2) 所有分支的下界是否满流
    # 使用第1种方法判定: 计算所有分支的下界之和
    # 用于后续判断是否可行流(满足下界,且下界之和等于最大流)
    S = np.sum([e['low'] for e in dg.es])
    ret = abs(S - max_flow) < 1.0e-2
    # 有上下界的最大流分配成功!!!
    if ret:
        # 得到计算的流量值
        for e in dg.es:
            e['flow'] = flow_map[e.index]
        # 记录拆分的分支
        edges = [e.index for e in dg.es if e['lu']>-1]
        # 记录已处理的分支
        aset = set()
        for i in edges:
            e = dg.es[dg.es[i]['lu']]
            if e.index not in aset:
                e['flow'] = e['low'] + flow_map[e.index]
                aset.add(e.index)
    # 删除属性数据
    del dg.es['capacity']
    del dg.es['lu']
    # 删除超级源汇点以及虚拟节点s0
    dg.delete_vertices([s0, ss, tt])
    return ret
# 最大流的固定风量分配算法
def fixed_max_flow(dg, Q, fixTotalQ=False):
    sn = add_virtual_source(dg)
    tn = add_virtual_target(dg)
    if sn < 0 or tn < 0:
        return False
    else:
        for i in dg.incident(sn, mode=OUT):
            dg.es[i]['weight'] = 0.0
        for i in dg.incident(tn, mode=IN):
            dg.es[i]['weight'] = 0.0
            
        #增加3个属性数据
        dg.es['low'] = 0.001
        dg.es['high'] = DBL_MAX #表示无穷大
        dg.es['flow'] = 0.0
        #分配固定风量
        for e in dg.es:
            if e['weight'] > 0:
                e['low'] = e['weight']
                e['high'] = e['weight']
        #调用有上下界的最大流算法
        ret = bound_max_flow(dg, sn, tn, Q, fixTotalQ)
        if ret:
            # 固定分量可分配,则分配所有分支的风量
            for e in dg.es:
                e['weight'] = e['flow']
        #删除属性数据
        del dg.es['low']
        del dg.es['high']
        del dg.es['flow']
        #删除虚拟源点和汇点,以及它们的关联分支
        dg.delete_vertices([sn, tn])
        return ret
#基于最大流的固定风量分配算法
def ifq(vnet):
    dg = vnet.graph()
    # 增加权重属性数据
    dg.es['weight']=0.0
    # 构造分支权重
    build_fixq_weight(vnet)
    # 调用基于最大流的固定风量分配算法
    ret = fixed_max_flow(dg, dg['Q'], False)
    if ret:
        # 分配计算的风量
        for e in dg.es:
            e['q'] = e['weight']
    # 删除权重属性
    del dg.es['weight']
    return ret
#构造边的权重--用于找最小生成树
def build_tree_weight(vnet, speedUp):
    # 计算分支权重
    dg = vnet.graph()
    # 放大固定风量分支和风机分支(作为余支)
    # r+1只是一个小技巧,避免0乘的情况
    # 依次给图中的其他分支分配权重
    for e in dg.es:
        r, q = R(e), e['q'] # 总风阻(包括调节风阻)
        c = (1.0+r)*q if speedUp else r+1.0
        # 风机分支
        if is_fan_edge(e):
            c = c*LARGE_COEFF
        # 固定风量分支
        elif is_fix_edge(e):
            c = c*LARGE_COEFF*10
        # 虚拟大气分支
        elif e.index == vnet.vAir():
            c = c*LARGE_COEFF*100
        e['weight'] = c
#找最小生成树
def find_mst(vnet, speedUp):
    dg = vnet.graph()
    #增加权重属性数据
    dg.es['weight']=0.0
    #构造分支权重
    build_tree_weight(vnet, speedUp)
    #igraph的最小生成树算法需要传递一个数组,用来表示分支的权重
    weight = [e['weight'] for e in dg.es]
    tree = dg.spanning_tree(weights=weight, return_tree=False)
    #删除权重属性
    del dg.es['weight']
    #返回找到的树边
    return tree
#通过pred属性得到路径
def node_path_from_pred(dg, s, t, pred_name='pred'):
    #记录路径(用节点表示)
    pred = dg.vs[pred_name]
    path = [t]
    u = t
    while pred[u] >= 0 and u != pred[u]:
        u = pred[u]
        path.append(u)
    if len(path) > 1:
        if path[-1] != s:
            path = []
        else:
            path.reverse()
    else:
        path = []
    return path
# 节点路径转换成分支路径
def node_path_to_edge_path(dg, node_path, directed=False):
    if len(node_path) < 2:
        return []
    # print node_path
    edge_path = []
    for u,v in zip(node_path[:-1], node_path[1:]):
        i = dg.get_eid(u, v, directed=directed, error=False)
        if i < 0:
            continue
        else:
            edge_path.append(i)
    return edge_path
def edge_path_from_pred(dg, s, t, directed=True, pred_name='pred'):
    # 通过先驱节点返回节点路径
    path = node_path_from_pred(dg, s, t, pred_name=pred_name)
    # 节点路径转换成分支路径
    return node_path_to_edge_path(dg, path, True)

# 深度优先搜索
def dfs_path(g, s, t):
    #igraph的python接口没有实现dfs,暂时用bfs替代
    ns, layer, pred = g.bfs(s, mode=OUT)
    g.vs['pred'] = -1
    for v in g.vs:
        v['pred'] = pred[v.index]
    # print ns, layer, pred
    # 返回节点路径
    return node_path_from_pred(g, s, t)
#找独立回路
def find_cycle(vnet, speedUp):
    # 查找最小生成树
    tree = find_mst(vnet, speedUp)
    if len(tree) == 0:
        return []
    dg = vnet.graph()
    
    # 通过一些分支构造子图
    # subgraph_edges
    # 通过一些节点构造子图
    # induced_subgraph
    # 构造一棵树(树是一颗子图)
    # 查找所有的余支(通过集合的"差集"运算得到余支集合)
    edges_set = set([e.index for e in dg.es])  #所有分支集合
    tree_set = set(tree)  #树支集合
    left_tree_set = edges_set.difference(tree_set) #余支集合
    # 将树边转换成无向图
    g = dg.subgraph_edges(list(tree))
    g.to_undirected()
    # print g
    #独立回路数组
    cl = []
    # 对每个余支搜索回路
    for i in left_tree_set:
        # 无向图dfs搜索,v->u的路径(只有1条路径)
        u, v = dg.es[i].tuple
        node_path = dfs_path(g, v, u)
        # print v, '-->', u, ': ', node_path
        #节点路径转换为分支路径
        path = node_path_to_edge_path(dg, node_path, False)
        if len(path) == 0:
            continue
        else:
            # 构造独立回路对象(以余支i作为回路的基准分支)
            cycle = IndependentCycle(dg, i, path)
            # 添加回路到链表中
            cl.append(cycle)
    return cl
#判断数据的是否满足要求精度
def is_meet_error_precise(dvalues, precise):
    if len(dvalues) == 0:
        return True
    else:
        return np.max(np.abs(dvalues)) <= precise

# 网络解算迭代过程.
# 基本过程: (1)查找回路 (2)回路迭代
# @param vnet 通风网络.
# @param count 迭代次数.
# @param speedUp 是否进行加速.
# @note "加速"顾名思义是为了加快收敛过程,查找最小生成树的时候利用|r*q|代替r作为分支的权重
# @note 数学上已证明该技巧可以显著的加快迭代收敛,参考<<流体网络理论>>.
def iterate(vnet, count, speedUp, bCrackFan = True):
    dg = vnet.graph()
    # 根据最小生成树查找独立回路
    cl = find_cycle(vnet, speedUp)
    # for cycle in cl:
    #   e = cycle.getBaseEdge()
    #   print 'e%s\t' % (e['id'])
    # print
    if len(cl) == 0:
        return False
    else:
        # 迭代计算
        ret = False
        k = 0
        while not ret and k <= count:
            DQ, DH = [], []
            for i,cycle in enumerate(cl):
                # 每个回路进行迭代
                each_ret, dq, dh = cycle.iterate(bCrackFan)
                # 回路不参与迭代!!!
                if not each_ret:
                    # e = cycle.getBaseEdge()
                    # print '第', k, '次迭代回路', e['id'], '不参与迭代计算'
                    continue
                # e = cycle.getBaseEdge()
                # print('第%d个回路(基准分支:e%s)第%d次迭代: dq=%.3f, df=%.3f' % (i+1, e['id'], k+1, dq, dh))
                DQ.append(dq)
                DH.append(dh)
            # 判断精度是否满足要求(风量精度和阻力精度)
            if is_meet_error_precise(DQ, dg['q_precise']) and is_meet_error_precise(DH, dg['h_precise']):
                ret = True
                break
            else:
                # Q, H = vnet.totalFlow(), vnet.totalH()
                # print k,Q,H, np.max(np.abs(DQ)), np.max(np.abs(DH))
                # print k,Q,H, DQ, DH
                k = k + 1
        # print '实际的迭代次数:',k-1
        # 销毁回路对象的内存
        del cl
        return ret

'''
先迭代20次，再迭代maxCount-20次,总迭代次数maxCount
'''
def __vno1(vnet):
    dg = vnet.graph()
    # 迭代计算结果
    ret = False
    if dg['maxCount'] < 20:
        # 按给定次数迭代
        ret = iterate(vnet, dg['maxCount'], False, True)
    else:
        # 先迭代20次(不加速speedUp=False)
        ret = iterate(vnet, 20, False, True)
        # 如果20次后还不收敛(迭代加速speedUp=True)
        if not ret:
            ret = iterate(vnet, dg['maxCount']-20, True, True)
    # 返回迭代结果
    return ret

'''
按照20次一个周期进行迭代.
'''
def __vno2(vnet):
    dg = vnet.graph()
    # 总的迭代次数
    totalCount = 0
    # 迭代结果
    ret = True
    while True:
        # 是否加速(如果20次计算后仍不收敛,开始加速)
        speedUp = (ret == False)
        # 迭代解算
        ret = iterate(vnet, 20, speedUp, True)
        # internal_print_network(vnet, '第%d次迭代风量' % k)
        # 总的迭代次数累加
        totalCount = totalCount + 20
        if ret or totalCount > dg['maxCount']:
            break
    # 返回迭代结果
    return ret

#网络解算(需要用户自己去添加虚拟源汇点)
def vno(vnet):
    dg = vnet.graph()
    # 初始化风量(包括固定风量),利用最大流固定风量分配算法分配
    # 传统的风网解算程序由于没有解决"风机特性曲线振荡", "迭代加速"问题,要求风量分配要尽可能接近真实风量
    # 本算法无此限制条件,风量分配的结果只要保证风量守恒即可!
    # internal_print_network(vnet, '初始风量(Q=%.1f)' % dg['Q'])
    if not ifq(vnet):
        print('初始化固定风量失败')
        return False
    else:
        # 添加虚拟大气分支
        vnet.addAirEdge()
        # 迭代计算(先迭代20次,然后迭代加速maxCount-20次)
        # ret = __vno1(vnet)
        # 迭代计算(按20次一个周期进行迭代)
        ret = __vno2(vnet)
        if ret:
            vnet.caclNodePressures()
        # 删除虚拟大气分支
        vnet.delAirEdge()
        # internal_print_network(vnet, '初始风量(Q=%.1f)' % dg['Q'])
        # 返回迭代结果
        return ret
# 通风网络解算(内部自动添加虚拟源汇)
def vno2(vnet):
    vnet.addVirtualST()
    # print '虚拟源点:',vnet.vSource(),  '虚拟汇点:',vnet.vTarget()
    # 通风网络解算
    ret = vno(vnet)
    # 删除虚拟源汇
    vnet.delVirtualST()
    return ret

# 回路阻力平衡参数
def cycle_pressures(vnet, debug=False):
    dg = vnet.graph()    
    # 添加虚拟大气分支
    vnet.addAirEdge()
    # 根据最小生成树查找独立回路
    cl = find_cycle(vnet, False)
    # 循环计算回路阻力代数和
    DH = []
    f = open('cycles.txt', 'w')
    f.write('编号\t余枝\t树枝\n')
    for i, cycle in enumerate(cl):
        dh = cycle.pressureDrop()
        if debug:
            baseEdge = cycle.getBaseEdge()
            path = ['e%s' % (e['id']) for e in cycle.getPath() if not is_zero_edge(e)]
            dr = -dh/(baseEdge['q']*baseEdge['q'])
            print('回路C%d={e%s,%s}\tΔP%d=%.4f\tΔr=%.6f(基准分支)' % (i+1, baseEdge['id'], ','.join(path), i+1, dh, dr))
            f.write('%d\te%s\t%s\n' % (i+1, baseEdge['id'], ' '.join(path)))
        DH.append(dh)
    # 删除虚拟大气分支
    vnet.delAirEdge()
    f.close()
    return DH

def build_pressure_weight(vnet):
    dg = vnet.graph()
    #计算分支权重
    # 以分支的阻力作为权重
    for e in dg.es:
        e['weight'] = f0(e)
def spfa(dg, s):
    # 最短路径搜索
    # print s
    dists = dg.shortest_paths_dijkstra([s], weights='weight')
    if len(dists) == 0 or len(dists[0]) == 0:
        return
    dists = dists[0]
    # 遍历所有节点
    for v in dg.vs:
        if not np.isinf(dists[v.index]):
            v['dist'] = dists[v.index]
def back_spfa(dg, s):
    # 构造反向图(将原图分支反向)
    # print [(e.source, e.target, int(e['weight'])) for e in dg.es]
    # fix by dlj 2017/1/19: 分支增加_do_mrp属性,用于屏蔽大气分支,避免找最大阻力路线出现单向回路
    # 只选择_do_mrp属性为True的分支(排除大气分支)
    seq = dg.es.select(do_mrp_eq=True)
    edges = [(e.target, e.source) for e in seq]
    weight = [e['weight'] for e in seq]
    g = Graph(directed=True)
    g.add_vertices(len(dg.vs))
    g.add_edges(edges)
    g.es['weight'] = weight
    g.vs['dist'] = DBL_MAX
    # 执行spfa算法
    spfa(g, s)
    # 设置节点的dist属性
    dg.vs['dist'] = g.vs['dist'][:]

# BFS和Dijkstra算法都是A*算法的一种(f(n)=g(n)+h(n))
# 其中BFS的h(n)=0,没有利用启发式信息,盲目搜索,而Dijkstra也是如此
# A*算法和Dijkstra算法都属于一类,原则上不能应用于有负权重分支的图
# 但本程序假设除了虚拟分支以外所有的权重都为负值,所以也是可以适用的
def a_star(dg, s, t, k=1):
    # 用于priority_queue的比较函数.
    # [知识点]:K最短路径
    # 节点的估价函数值,决定A*算法的搜索方向
    # 对于A*算法,估计函数=当前值+当前位置到终点的距离成本
    # 即f(p)=g(p)+h(p),其中g(n)表示
    # 对于k最短路径算法来说,g(p)为当前从到p所走的路径的长度
    # h(p)为点p到t的最短路的长度
    class Node:
        def __init__(self, u, pred, g=0, h=0):
            self.u = u
            self.pred = pred
            self.g = g
            self.h = h
        def __lt__(self, other):
            f1 = self.g + self.h
            f2 = other.g + other.h
            # print abs(f1-f2), self.g, other.g, f1, f2
            if abs(f1-f2) < 1.0e-2:
                # print self.g, '<', other.g
                return self.g < other.g
            else:
                # print f1, '<', f2
                return f1 < f2
    if s == t:
        return False
    pq = PriorityQueue()
    cnt = 0
    dg.vs[s]['pred'] = -1
    node = Node(s, -1, 0, dg.vs[s]['dist'])
    pq.put(node)
    while not pq.empty():
        node = pq.queue[0]
        u = node.u
        dg.vs[u]['pred'] = node.pred
        # print 'pop: u=%d f=%.3f g=%.3f  h=%.3f' % (node.u, node.g+node.h, node.g, node.h)
        pq.get() # 类似pq.pop()
        if u == t:
            cnt = cnt + 1
        if cnt == k:
            break
        for i in dg.incident(u, mode=OUT):
            e = dg.es[i]
            v = e.target
            dg.vs[v]['pred'] = u
            child_node = Node(v, u, node.g+e['weight'], dg.vs[v]['dist'])
            # print '\tpush: v=%d f=%.3f g=%.3f h=%.3f' % (child_node.u, child_node.g+child_node.h, child_node.g, child_node.h)
            pq.put(child_node)
    return cnt == k
# 权重逆转(取负值 或 用一个较大的数减去权重值)
def reverse_weight(dg):
    max_weight = 0.0
    for e in dg.es:
        e['weight'] = max_weight - e['weight']

# 搜索2个节点之间的第k最大阻力路线
def mrp_k(vnet, s=-1, t=-1, k=1, bReverse=True):
    if k<=0:k=1
    dg = vnet.graph()
    # fix by dlj 2017/1/19: 分支增加_do_mrp属性,用于屏蔽大气分支,避免找最大阻力路线出现单向回路
    dg.es['do_mrp'] = True
    if vnet.vAir() > -1:
        dg.es[vnet.vAir()]['do_mrp'] = False

    # 如果始节点的值"无效"(INVALID),那么默认始节点为虚拟的源点
    if s < 0:
        s = vnet.vSource()
    # 如果末节点的值"无效"(INVALID),那么默认末节点为虚拟的汇点
    if t < 0:
        t = vnet.vTarget()
    # 如果k小于等于0,则默认k=1(也就是搜索第1最大阻力路线)
    if k <= 0:
        k = 1
    # 节点增加dist属性和pred属性
    dg.vs['dist'] = DBL_MAX
    # 节点的前驱节点(父节点),用于后续计算路径
    dg.vs['pred'] = -1
    # 分支增加weight属性
    dg.es['weight'] = 1.0
    
    # 以分支阻力为权重
    build_pressure_weight(vnet)
    # 权重取负值,将求最大权重相关的问题(如最长路径,最大阻力路线等)转换成最短路径问题)
    if bReverse:
        reverse_weight(dg)
    # 计算节点到末节点t的最短路径
    back_spfa(dg, t)
    # 使用A*算法搜索k最短路径
    k_shortest_path = []
    if a_star(dg, s, t, k):
        k_shortest_path = node_path_from_pred(dg, s, t)
    # 节点路径转换成分支路径
    k_shortest_path = node_path_to_edge_path(dg, k_shortest_path, True)
    # 过滤虚拟分支
    k_shortest_path = vnet.filterVirutalEdges(k_shortest_path)
    
    # 删除weight属性
    del dg.es['weight']
    # 删除dist属性
    del dg.vs['dist']
    # 删除pred属性
    del dg.vs['pred']
    # fix by dlj 2017/1/19: 分支增加_do_mrp属性,用于屏蔽大气分支,避免找最大阻力路线出现单向回路  
    # 删除_do_mrp属性
    del dg.es['do_mrp']
    # 返回k最大阻力路线
    return k_shortest_path

# 搜索两个节点之间的最大阻力路线
def mrp(vnet, s=-1, t=-1, bReverse=True):
    return mrp_k(vnet, s, t, 1, bReverse)

# 搜索整个通风系统的n条最大阻力路线
def mrp_n(vnet, n=1, bReverse=True):
    if n<=0:n=1
    paths = []
    for k in range(1, n+1):
        # 搜索第k最大阻力路线
        path = mrp_k(vnet, k, bReverse)
        if len(path) == 0:
            break
        else:
            paths.append(path)
    return paths

# 搜索单个风机的第K最大阻力路线
def mrp_fan_k(vnet, fan, k=1, bReverse=True):
    if k<=0:k=1
    dg = vnet.graph()
    if fan < 0 or fan >= len(dg.es):
        return []

    s, t = vnet.vSource(), dg.es[fan].source
    if t < 0:
        return []
    else:
        # 搜索最大阻力路线
        maxP = mrp_k(vnet, s, t, k, bReverse)
        if len(maxP) > 0:
            maxP.append(fan)
            return maxP
        else:
            return []

# 搜索单个风机的最大阻力路线
def mrp_fan(vnet, fan, bReverse=True):
    return mrp_fan_k(vnet, fan, 1, bReverse)

# 不考虑压差限制,搜索单个风机的n条最大阻力路线
def mrp_fan_n(vnet, fan, n=1, bReverse=True):
    if n<=0:n=1
    dg = vnet.graph()
    if fan < 0 or fan >= len(dg.es):
        return []

    paths = []
    for k in range(1, n+1):
        # 搜索第k最大阻力路线
        path = mrp_fan_k(vnet, fan, k, bReverse)
        if len(path) == 0:
            break
        else:
            paths.append(path)
    return paths

# 给定压差h限制,搜索单个风机的n条最大阻力路线
def mrp_fan_h(vnet, fan, n=1, minDeltaH=0, maxDeltaH=DBL_MAX):
    # 调整负风量分支的风流方向
    # vnet.adjustNegativeEdges()
    if n<=0:n=1
    dg = vnet.graph()
    if fan < 0 or fan >= len(dg.es):
        return []

    s, t = vnet.vSource(), dg.es[fan].source
    if t < 0:return []
    # 先搜索单个风机的最大阻力路线
    maxP = mrp(vnet, s, t)
    if len(maxP) == 0:
        return []

    # 计算最大阻力路线的总阻力
    maxH = vnet.hPath(maxP)
    # 记录所有的最大阻力路线
    paths = []
    # 搜索成功一次就累加一次,直到最大阻力路线的条数达到n个
    k = 1
    while len(paths) < n:
        path = mrp_k(vnet, s, t, k)
        if len(path) == 0:
            break
        else:
            # 搜索成功一次就累加一次
            k = k+1            
            # 计算次最大阻力路线的总阻力
            H = vnet.hPath(path)
            # 只考虑总阻力的差值大于tol的次最大阻力路线
            deltaH = abs(maxH - H)
            # print 'path=',path, 'δH=',deltaH
            if deltaH >= minDeltaH and deltaH <= maxDeltaH:
                # 搜索成功后将风机分支添加路径中
                path.append(fan)                
                paths.append(path)
    return paths

# 搜索所有风机的最大阻力路线(没有太大的实际用途!)
def mrp_fans(vnet, bReverse=True):
    PP = []
    fan_edges = vnet.fanEdges()
    # 查找每台风机的最大阻力路线,并打印出来
    for fan in fan_edges:
        # 查找到风机的最大阻力路线
        P = mrp_fan_k(vnet, fan, bReverse)
        PP.append((fan, vnet.filterVirutalEdges(P)))
    return PP

# 查找多个回风井(汇点)的公共分支
def apm_common_edges(vnet):
    dg = vnet.graph()
    #查找所有出度等于0的节点
    nodes = dg.vs.select(_outdegree=0)
    if len(nodes) == 0:
        return []
    elif len(nodes) == 1 and nodes[0] == vnet.vTarget():
        nodes = [dg.es[i].source for i in dg.incident(nodes[0], mode=IN)]
    if len(nodes) < 2:return []
    # 所有通路共有的公共分支
    common_edges = []
    # 多个节点依次进行反向深度优先搜索
    arc_set = set()
    for s in nodes:
        ns, layer, pred = g.bfs(s, mode=IN)
        for u in ns:
            if pred[u] < 0:continue
            i = dg.get_eid(u, pred[v], directed=True, error=False)
            if i < 0:
                continue
            if i not in arc_set:
                arc_set.add(i)
            else:
                common_edges.append(i)
    return common_edges

'''
采用通路法对一台风机进行调节(Adjust Path Method).
该函数只是一个样板函数,没有完全实现,用来大概描述通路法的思想,仅作参考(代码存在错误)!!!
@param  vnet  [in] 单一源汇通风网络(连通且无单向回路).
@param  fan   [in] 一台风机所在的分支
@note 通路法调节的主要思路(来自刘老师的答疑):\n
      "以风机为单位,一个风机一个风机的调.\n
      注意原则:最大阻力路线不能调,次最大调完后,次最大上的道也不能再调.\n
      总之,后调的不能把前面的最大阻力路线给变了!".
@note 我的补充思路:\n
      针对一条风机进行调节,搜索到最大阻力路线后,计算该条路线的总阻力maxH.\n
      将路线上的所有分支都着色(过滤掉),然后再搜索最大阻力路线(也就是次最大阻力路线),\n
      计算次最大阻力路线的总阻力,与maxH的差值就是要调节的量,然后根据每条分支的"可调属性",\n
      查找路线上的所有可调分支,然后让用户选择"调哪几条分支,调多少"(这是一个交互过程).\n
      用户调节完之后返回,再将"次最大阻力路线"上的分支着色过滤,循环上述步骤,直到不能再调为止.
@note 在调用函数前,必须将网络转换成单一源汇网络!!!
'''
def apm_sample(vnet, fan):
    # 预处理检查
    dg = vnet.graph()
    if fan < 0 or fan >= len(dg.es):
        return False
    s, t = vnet.vSource(), dg.es[fan].source
    if t < 0:
        return False

    # 首先搜索单个风机的最大阻力路线
    maxP = mrp(vnet, s, t)
    if len(maxP) == 0:
        return False

    # 计算最大阻力路线的总阻力
    maxH = vnet.hPath(maxP)
    # 压差限制
    minDeltaH, maxDeltaH = 0.1, DBL_MAX
    # 开始进行调节
    while True:
        # 搜索压差大于minDeltaH的最大阻力路线
        paths = mrp_fan_h(vnet, fan, minDeltaH, maxDeltaH, n=1)
        if len(paths) == 0:
            break
        # 计算当前最大阻力路线的阻力
        P = paths[0]
        H = vnet.hPath(P)

def print_node_pressure(vnet, msg='通风网络节点压力计算结果', nodes=[]):
    dg = vnet.graph()
    print('----------%s----------' % msg)
    V = dg.vs
    if len(nodes) > 0:
        V = nodes
    for i, v in enumerate(V):
        print('p(v%s)=%.3f\t' % (v['id'], v['p']))
        if (i+1)%3 == 0:
            print()
    print()    

def bool2str(x):
    if x:
        return '是'
    else:
        return '否'
    
def print_network(vnet, msg='通风网络分支风量分配结果', edges=[], debug=False):
    dg = vnet.graph()
    if debug:
        print('---------- igraph内部详细信息 ---------')
        print(dg)
    print('---------- %s ---------' % msg)
    E = dg.es
    if len(edges) > 0:
        E = edges
    print('序号\t分支\t始节点\t末节点\t风阻\t风量\t压力')
    for i, e in enumerate(E):
        u, v = dg.vs[e.source], dg.vs[e.target]
        print('%d\te%s\tv%s\tv%s\t%.4f\t%.4f\t%.4f' % (i+1,e['id'], u['id'], v['id'], e['r'], e['q'], f0(e)))
        # print 'e%d=(v%d,v%d)\tr(e%d)=%.3f\tq(e%d)=%.3f\th(e%d)=%.3f' % (e.index,u.index, v.index, e.index,e['r'],e.index,e['q'],e.index,f0(e))
    print('总分支个数:%d\t总节点个数:%d\t总风量:%.2f\t总阻力:%.2f\t总风阻:%.4f' % (len(dg.es), len(dg.vs), vnet.totalFlow(), vnet.totalH(), vnet.totalR()))
    e1, e2 = vnet.minMaxQ()
    print('最小风量:q%s=%.4f\t最大风量:q%s=%.4f' % (e1['id'], e1['q'], e2['id'], e2['q']))
    # 打印固定风量分支
    print('固定风量分支:%s' % (','.join(['e%s' % (dg.es[i]['id']) for i in vnet.fixQEdges()])))
    print('风机分支:%s' % (','.join(['e%s' % (dg.es[i]['id']) for i in vnet.fanEdges()])))
    # 判断回路时候平衡(并打印内部的回路)
    ret = np.max(np.abs(cycle_pressures(vnet, debug=True))) < dg['h_precise']
    # print(np.max(np.abs(cycle_pressures(vnet, debug=True))),dg['h_precise'])
    print('回路是否平衡:%s' % bool2str(ret))
    print('')

def print_mrp_fans(vnet):
    fan_edges = vnet.fanEdges()
    # 查找每台风机的最大阻力路线,并打印出来
    for fan in fan_edges:
        # 查找到风机的最大阻力路线
        P = mrp_fan_k(vnet, fan)
        dg = vnet.graph()
        print('最大阻力路线P(e%s)={%s}' % (dg.es[fan]['id'], ','.join(['e%s' % (dg.es[i]['id']) for i in P if not is_zero_edge(dg.es[i])])), 'h=%.2f' % (vnet.hPath(P)))

# 打印一条最大阻力路线
def print_mrp_fan(vnet, P):
    dg = vnet.graph()
    print('P={%s}' % (','.join(['e%s' % (dg.es[i]['id']) for i in P if not is_zero_edge(dg.es[i])])), 'H=%.2f' % (vnet.hPath(P)))

# 打印最大阻力路线
def print_mrp(vnet):
    print_mrp_fan(vnet, mrp(vnet))